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We present new molecular dynamics results for the pressure of the pure hard disk fluid up to the hexatic transition (about reduced density 
0.9). The data combined with the known virial coefficients (up to -Bio) are used to build an equation of state, to estimate higher-order 
virial coefficients, and also to obtain a better value of Bio ■ Finite size effects are discussed in detail. The "van der Waals-like" loop reported 
in literature in the vicinity of the fluid/hcxatic transition is explained by suppressed density fluctuations in the canonical ensemble. The 
inflection point on the pressure-density dependence is predicted by the equation of state even if the hexatic phase simulation data are 
not considered. 
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1 Introduction 

The hard disk fluid is the simplest model of various surface systems (surfactants, adsorption of molecules 
on smooth surfaces). It has been intensively studied by a number of methods, of which we are interested 
in Monte Carlo and molecular dynamics simulations ^[^llSlll]) virial coefficients [SHH]) and equations of 
state (EOS) [TllcM ITOlllTMOICTIliau^ . Recently, peculiarities of phase transitions from fluid to 

hexatic to solid are frequently studied [THIH]. 

In this work we focus mainly on the fluid region. The fluid phase changes at density of about p c = 
0.8995(23) (p denotes the reduced number density, p = Na 2 /A, where N is the number of disks, a their 
diameter, and A the system area) to the hexatic phase [iTHl2*U| IT%] ■ 



2 Methods 

2.1 Simulation methodology 

2.1.1 Molecular dynamics code. To obtain accurate data on the EOS, we use standard molecular 
dynamics (MD) simulations in the microcanonical ensemble and tetragonal (square) periodic boundary 
conditions with zero total momentum (MD-NVE; we keep the usual abbreviation NVE for constant Number 
of particles- Volume-Energy even if the volume is here replaced by the area) . Our MD program combines 
the ideas of the linked-cell list method both in space and time and is highly optimized. Details on the 
code are given in our previous paper In the present study, we use N = 4000 and 9000, and at higher 
densities 16 000, 25 000 and 50 000 disks of unit diameter in a square periodic box. 

2.1.2 Compressibility factor. The compressibility factor Z = pA/(NkT) = p/(pkT) (p denotes pres- 
sure, k the Boltzmann constant, and T absolute temperature) can be calculated from MD simulations in 
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Figure 1. Weights w for mixing both formulas for Z in equation Jl) and the fit bounded by to < 2. 



two ways PJE] . One formula uses the virial of force 

r-i 

IT 2E kin (t 2 - tx) 



N-l 1 . 



where the sum is over all collisions occurring during time interval (ti , ^2) ? ^ v ij is the change of velocity 
of both colliding disks of mutual position r^, and E k i n is the kinetic energy (constant in the MD-NVE 
simulation). The alternative formula uses the collision rate 



Z Ta t xs (t 1 ,t2) = \ + MN)J Z — - — V 1, (2) 
rateV 1, 2) t jy ,y 2DNE kiD t 2 - h f- ' ' K J 

where YlteCti t 2 ) ^ ^ s ^ e numDer of collisions in time interval (ti,t 2 ) and |2j 

= r[(g(jy-i) + i)/2] 

71 ' T[D(N — l)/2](DN/2) l l 2 ' {) 

where D = 2 is the dimensionality. 

Both formulas give within statistical inaccuracies the same values, but their statistical errors differ. Our 
final result is therefore a weighed average of both formulas, 

Z = wZ viT + (1 - w)Z Tate , (4) 

where the mixing weight w = w{p) is a function of density. Methodology for determining the optimum 
function w has been explained in detail in |21j . The simulation results on w{p) are shown in figure Q along 
with a function fitted to the data and bounded by w < 2 to avoid too extrapolated values. This function 
was used in final analysis of MD data to avoid bias. 
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Table 1. Molecular dynamics results 
on the compressibility factor, cor- 
rected for the finite size effects; for 
p < 0.88 by equation ® 5 for p > 0.88 
by Z(lfN) linear extrapolation. 



p 


z 


cr(Z) 


0.40 


2.1514393 


0.0000052 


0.45 


2.4276680 


0.0000065 


0.50 


2.7601235 


0.0000079 


0.55 


3.1647878 


0.0000098 


0.60 


3.663691 


0.000012 


0.65 


4.287926 


0.000015 


0.70 


5.082362 


0.000019 


0.75 


6.113391 


0.000026 


0.80 


7.476491 


0.000036 


0.83 


8.494891 


0.000050 


0.84 


8.866011 


0.000059 


0.85 


9.245785 


0.000072 


0.86 


9.621609 


0.000082 


0.87 


9.96782 


0.00013 


0.88 a 


10.2309 


0.0003 


0.89 


10.3176 


0.0011 


0.90 b 


10.2059 


0.0011 


A 


compromise 


between 



10.23095(26) by and 
10.23055(54) by Z(l/N) ex- 
trapolation. 

b This value may be affected by fi- 
nite size effects. 

2.1.3 Start and equilibration. The initial configuration was a tetragonal crystal with random velocities 
assigned to particles. This setup creates a disorder within a few collisions. At higher densities, a locally 
hexagonal arrangement gradually develops. Therefore a period of equilibration follows until the compress- 
ibility factor reaches a constant value with fluctuations. At higher densities and for a large system this 
takes a long time. 



2.2 Finite size effects 

2.2.1 Fluid region. The finite-size errors are in our pseudoexperimental setup at low densities several 
times larger than the statistical errors. E.g., at the lowest density simulated, p = 0.4, one would need 
as many as N = 160 000 particles to guarantee the same systematic finite-size error as the statistical 
error of tabled It is more efficient to simulate a smaller system and to correct for the finite-size error to 
obtain the thermodynamic limit [N — > oo) because correlation times and therefore convergence in large 
two-dimensional systems is slow. At about p = 0.85 both errors become comparable for N = 20 000. In a 
vicinity of the hexatic transition the finite-size errors become large again and difficult to determine. 
The correction procedure can be divided into three steps [2*Tll2*2*] : 

MD-NVE — > NVT — > pVT — ► thermodynamic limit, 

where the first step is already included in formulas (TjQ) and P|). and ^VT stands for the grand-canonical 
ensemble. 

The largest source of inaccuracy for hard-body systems is the second step, NVT — ► pVT |2MU24| . This 
"ensemble correction" is caused by suppressed fluctuations of density (number of particles) in the canonical 
ensemble and is given by |2*T| l23] 

-l 

+ 0(N- 2 ), (5) 



ry H " f 

■"NVT = — T 



p dp 
2N ~dp 2 



dp 
dp 



where indices pVT and NVT refer to the respective expectation values in the grand-canonical and canonical 
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ensembles at the same number density (i.e., (iV^vT = -Wnvt)- This term requires knowledge of an EOS 
which is not known in advance. The correction procedure is therefore self-consistent (we start with no 
correction, fit the EOS, calculate the correction and better data, etc.), albeit rapidly converging. For 
p > 0.88, the accuracy of the EOS is not sufficient to calculate the second derivative needed in Q and we 
use linear extrapolation of Z in dependence on 1/N to zero. In a close vicinity of p c even this approach 
fails. 



2.2.2 Near the critical point. It was investigated in detail in JH] that the p(p) dependence has an 
inflection point at the fluid/hexatic transition, p = ZkTp = po — const x (p c — p) a as p — > p c , where 
a' = 4.55. Then the "ensemble error" (of opposite sign than the ensemble correction) becomes 

Znvt " Z ^ VT = J-^p 277 (6) 

and similarly for p > p c (likely with a different critical exponent a'). This term is positive for p < p c , 
negative for p > p c , and it diverges at p c . For p = 0.895 and N = 256 2 we obtain from Q the NVT 
pressure by 0.005 higher than the grand-canonical one which is in qualitative agreement with the value of 
0.01 by fl] (the statistical error is not provided, but it is likely less than 0.01). The correct value of the 
correction may be affected by higher order terms; apparently this must be true in a close vicinity of p c 
because the ensemble correction cannot diverge at finite N. Therefore the "van der Waals-like" loop [3] is 
an artifact caused by suppressed density fluctuations in the constant- N ensemble and it disappears if data 
are free of finite-size effects. 

Near the first-order phase transition (occuring e.g. for hard spheres) one observes a hysteresis caused 
by long-lived metastable states rather than a "loop"; a loop can arise at very small systems and for long 
runs so that an equilibrium is maintained. There is no true hysteresis near the continuous transition point, 
although of course the dynamics (and convergence) slows down. 



2.3 Equation of state 

2.3.1 Correlation of the data. The pseudoexperimental data were fitted to a polynomial in x = y/(l — 
y), where y = ^4hdP is the packing fraction and ^4hd is the disk area. The compressibility factor reads as 

i=o v y 7 

where Aq-A^ are determined so that virial coefficients B2-B5, known either analytically |26I27| or (S5 |28| ) 
with high precision, are exactly reproduced. Parameters Ai, 4 < i < k, are adjustable, but some of them 
may be zero. The number of degrees of freedom (number of all data minus number of adjustable parameters) 
is denoted nf ree . 

This choice needs some explanation. Terms y/(l — y) appear in several theories for both hard spheres and 
hard disks: the scaled particle theory [2H], Percus-Yevick and hypernetted-chain integral equations (3131) an d 
resummation of approximated virial coefficients (Carnahan-Starling equation) |31| . The main reason for 
this form is that it moderates the sharp increase of Z with increasing y as suggested in |32| (equation (13)) 
under name y-expansion (function y is here called x and should not be confused with packing fraction y). 

Another popular choice is a rational function which is closely related to the Pade approximant [HHj 
derived from a formal virial power series; the coefficients are determined by the least-square method. 
There is some experience with this method in our laboratory |34j . but detailed investigation including 
high densities found this method numerically unstable and less flexible; we have not checked this approach 
for hard disks, though. In addition, we have theoretical objections against interpretation of the pole (zero 
point of the denominator) of such a rational function as the random close packing of the hard sphere fluid 
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(frozen glass). This interpretation assumes that the EOS can be analytically continued beyond the freezing 
point which is not true 35,36,37 ; it can be continued with a limited precision only — the random close 
packing is an inaccurate concept in principle. With increasing precision of input data this inaccuracy may 
cause problems in fitting. In the hard disk system it is impossible to extrapolate beyond the continuous 
transitions (with inflection points on the pressure-density dependence) to get any random close packing 
and even this inaccurate reason for using rational functions becomes invalid. 

Another consequence of the nonanalyticity is that the radius of convergence of the virial series is less 
than or equal to the phase transition density (first-order freezing of the hard sphere fluid or continuous 
transition of the hard disk fluid). Consequently the approximately quadratic -B* dependence (for hard 
spheres) and linear dependence (for hard disks), related to expanded y/(l — y) terms, cannot extent to 
infinite n. It is in principle possible to determine the radius of convergence from the sequence of B n , but 
the available precision and the number of terms are not sufficient. 

With functional form (J7J), the standard objective function 



1 



n fr , 



"data 

E 

3=1 



+ 



10 

E 

j=6 



£>EOS _ ^MC 

3 a(Bf c ) 



(8) 



was minimized, where a stands for the standard error, Bj os is the virial coefficient calculated from the 
EOS and Bf c ± a(Bf c ) are virial MC data with standard errors 0E]. In other words, both the MD data 
on the compressibility factor and the MC data on the virial coefficient are correlated simultaneously. 

The value of s for an optimum fit is around unity provided that the input standard errors a are reliable, 
which is the case of our simulations where a is determined with accuracy (error of error) of a few per 
cent |21j . If s 3> 1 then the number of adjustable parameters is not sufficient do describe the data. On the 
other hand, one should not use more parameters than necessary because just noise would be fitted; the 
best test is to remove one parameter and to observe whether s significantly increases. 



2.3.2 Higher-order virial coefficients. Expanding equation (JJJ) in powers of density gives virial coef- 
ficients. Virial coefficients Bi for i < 5 are exactly reproduced, for 6 < % < 10 they are modified because 
their change is allowed by a simultaneous fit (within statistical errors), and for i > 10 they are predicted. 



3 Results 

3.1 Molecular dynamics 

First, we have checked that both the virial © and rate © routes to the compressibility factor are 
equivalent. The differences are within combined statistical errors the same and only for two systems with 
N = 9000 and p = 0.7 and p = 0.75 slightly exceed one standard deviation. In fact, because of correlations 
in the data, the difference is at low and large densities much less than the combined standard error. 

In order to assess the role of the periodic errors, we fitted the total correlation function h(r) = g{r) — 1 
at large separations to attenuated oscillations, h(r) = He[Aexp(—Br)/r], where A and B are complex 
constants and Re denotes the real part. The value of Re(U) describes the decay of correlations and its 
typical value at simulation square size L = A 1 / 2 , Re(^4) exp[— ~Rje(B)L]/L, is the estimate of the periodic 
error. We found that Re(-B) = 0.22 for p = 0.88. Consequently the periodic error is negligible for all used N 
and p < 0.88, and therefore formula © or linear Z(l/N) dependence is sufficient to account for finite size 
effects. For p = 0.89 we found Re(-B) = 0.11 and N = 4000 may have a small periodic error about 2 • 10~ 5 ; 
the datum was nevertheless discarded. For p = 0.90 ~ p c it holds Re(B) = 0.044 and even N = 16 000 is 
barely sufficient. 

For p < 0.88, the data for different iV were corrected by © and a weighed average was taken; a check 
was made that the corrected data do match within statistical errors. Correction term © is not applicable 
for p > 0.89 because it contains the second derivative of the EOS and therefore the correction term is 
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Figure 2. The data and proposed EOS. Dotted line: p ma x = 0.88, dashed line: p ma x = 0.89, solid line: p ma x = 0.90. 

large and not available with sufficient precision. Linear extrapolation of Z(l/N) was used instead; the 
final results thus cease precision. Point p = 0.88 is a borderline between applicability of both approaches: 
The ensemble correction is large and its accuracy may affect the results, but the extrapolation gives less 
accurate data. 

The final corrected MD data are collected in tabled The data in the "difficult" region close to the phase 
transition agree well with recent extensive Monte Carlo data (HE] with the exception of density p = 0.9 
closest to p c where the N = 1024 2 result Z = 10.212 (HJ is significantly larger than our Z = 10.206; in 
this case one cannot assume linearity of the Z(l/N) dependence and our value is probably affected by 
higher-order finite size errors. 



3.2 Equations of state 

We present three best versions differing by the maximum density /9 max , number of fitted parameters and 
the value of the objective function s (JSJ), see figured Note that x = y/(l — y), where y is the packing 
fraction. 

Pmax = 0.88, s = 0.724: 

Z = l + 2x + 1.12801775x 2 + 0.00181895291 x 3 - 0.0526134737 x 4 + 0.0504951668 x 5 - 
0.0325433846 x 6 + 0.0133946531 x 7 + 0.00174265604 x 8 - 0.00944632202 x 9 + 0.00851111768 x 10 - 
0.0035963525 x 11 + 0.000577345106 x 12 - 1.06399127 • 10" 7 x 19 

/W = 0.89, s = 0.966: 

Z = l + 2x + 1.12801775x 2 +0.00181895291 x 3 - 0.0526134737 x 4 + 0.0504963915 x 5 - 
0.0325578581 x 6 + 0.0134816028 x 7 + 0.00129187484 x 8 - 0.00808881628 x 9 + 0.00669011963 x 10 - 
0.00250795961 x 11 + 0.000336036442 x 12 - 5.15282664 • 10~ 9 x 22 + 5.57730095 • 10~ 23 x 57 

Pmax = 0.90, a = 0.927; region p G [0.89, 0.90] of this equation may be affected by finite size effects: 

Z = l + 2x + 1.12801775x 2 +0.00181895291 x 3 - 0.0526134737 x 4 + 0.0504960168 x 5 - 
0.0325537792 x 6 + 0.0134578632 x 7 + 0.00140888182 x 8 - 0.00834273601 x 9 + 0.00694127367 x 10 - 
0.00262254723 x 11 + 0.000355746352 x 12 - 5.24672938 • 10" 9 x 22 + 5.88054639 • 10~ 23 x 57 



2008 14:20 Molecular Physics hdeos 



Hard disk fluid 7 




10 11 12 13 14 15 16 17 18 19 20 

n 

Figure 3. Higher-order virial coefficients as predicted by different EOSs. The best equations from Section lij, 21 are marked by * in the 
table and by solid lines, otherwide the dash size is proportional to the value of s (the longer dash, the worse fit). The key determines 
the equation: a digit (A=10, etc.) stands for one adjustable parameter (in addition to five parameters corresponding to Be to Bio), its 
value is the power of y with respect to the previous term. Error bars are standard errors propagated from standard errors of the MD 

data and known virial coefficients. 

It is interesting that both equations with p max > 0.89 predict to some extent the loop (with the "classical" 
critical exponent a' = 3) at the critical (fluid/hexatic) point, even if this is not the aim of the present 
work which focuses rather on the low-density region. Any extrapolation to p > p c should be done with 
caution because function p(p) is likely nonanalytical at p c . 

3.3 Virial coefficients 

The virial coefficients were determined from a number of EOSs, see figure The error of a predicted virial 
coefficient consists of a statistical error which is calculated from standard errors of input data and from a 
systematic (method) error. The former error ranges from typically 0.0014 (max. 0.002) for B* to typically 
0.02-0.05 (max. 0.07) for to typically 0.04-0.2 (max. 0.4) for £>* 6 . The latter error is apparently 
difficult to determine. We believe that the range of different predictions, also those with s > 1, gives a 
certain measure the systematic error; it also approximately matches the maximum statistical errors found 
in the used set od EOSs. We are rather pessimistic in determining this error; most EOSs with s < 1 (thin 
lines in figure |2J give several times less scattering in the virials predicted and only moderately increased 
statistical errors. All virial coefficients are collected in tabic [2j 

To verify the procedure, we repeated the calculations with Bio removed from the second sum of ©• The 
predicted value was B\q = 10.2210 ± 0.002, which is in agreement with the direct MC datum [H] (and in 
fact more accurate). The "best" value based on all available data (incl. B\q of f> ) is only slightly smaller, 
Bio = 10.2203 ± 0.002. In contrast, lower-order virials are accurate enough and including the MD data in 
the fit does not improve precision. 

4 Concluding remarks 

The proposed equations of state in the fluid region combine all available information — virial coefficients 
and simulation compressibility data. The equations may serve in perturbation theories. They are not meant 
as a replacement of physically-based (but less accurate) equations which are able to described more phases. 
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Tabic 2. Virial coefficients reduced by 
the disk area (packing fraction expan- 
sion), B? = BiA^j. The B G to B 9 val- 
ues are combined data of islisl - for Bio 
see ref. 0, Bn to B 16 arc our EOS-bascd 
predictions with conservatively estimated 
error bounds (not standard errors). 



n 


K 




2 


2 





3 


3.1280177516 


(rounded) 


1 


4.2578544562 


(rounded) 


5 


5.33689664 


0.00000064 


6 


6.3630259 


0.0000109 


7 


7.352077 


0.000028 


8 


8.318677 


0.000061 


9 


9.27234 


0.00027 


10 


10.2161 


0.0041 


10 a 


10.2203 


0.002 


10 b 


10.2210 


0.002 


11 


11.172 


0.010 


12 


12.132 


0.03 


13 


13.097 


0.06 


14 


14.053 


0.08 


15 


14.94 


0.21 


16 


15.7 


0.4 



a EOS-based prediction incl. Bio °f HI; 
the recommended value. 



b EOS-based prediction (-Bio of 5 not 
used) . 

The equations enable prediction of higher-order virial coefficients with no additional assumption on 
their order-dependence. The results also witness about the "law of complexity conservation": The value 
of the tenth virial coefficient can be obtained with comparable precision both directly by diagrammatic 
techniques [5] and by simulations. 

In order to obtain highly accurate MD data, it was necessary to take into account finite-size effects, which 
is especially peculiar close to the critical fluid/hexatic point. The "van der Waals-like" loop reported by 
several authors in this region can be semiquantitatively predicted by the concept of suppressed density 
fluctuations in the canonical ensemble. 
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